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Abstract 

We present valence electron Compton profiles calculated within the density-func- 
tional theory using the all-electron full-potential projector augmented-wave method 
(PAW). Our results for covalent (Si), metallic (Li, Al) and hydrogen-bonded ((H20)2) 
systems agree well with experiments and computational results obtained with other 
band-structure and basis set schemes. The PAW basis set describes the high-momen- 
tum Fourier components of the valence wave functions accurately when compared 
with other basis set schemes and previous all-electron calculations. 
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1 Introduction 



Compton scattering spectroscopy [1] is a method for obtaining direct infor- 
mation on electronic states in materials. From the scattering cross section one 
can extract the Compton profile J{pz), which, within the impulse approxi- 
mation [2], is a one-dimensional projection of the three-dimensional electron 
momentum density p{p). Compton scattering experiments are often analyzed 
with the help of electronic structure calculations, which are usually based on 
the framework of the band theory in the local-density approximation (LDA) of 
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the density-functional theory [3] (DFT). A variety of different methods can be 
used to represent single-particle states, charge density and potential in prac- 
tical calculations. The approximations made can have a large impact on the 
quality of the Compton profiles obtained. 

The projector augmented-wave method (PAW) introduced by P. E. Blochl [4] 
in 1994 is an elegant all-electron method generalizing both the pseudopoten- 
tial method and the linearized augmented-plane-wave method (LAPW). The 
method includes a well-defined linear transformation from soft pseudo-valence 
(PS) wave functions used in the actual calculation to all-electron (AE) valence 
wave functions. This is very practical for the calculation of electron momen- 
tum densities, because the AE wave functions with the rapid oscillations can 
thus be easily constructed. One then obtains properly also the high-momentum 
Fourier components of the wave functions. The PAW method has already been 
applied to many computational problems that require knowledge of the AE 
wave functions or the full charge density near the nuclei, for example, to the 
calculation of magnetic hyperfine parameters [5], electric-field gradients [6], 
and optical properties of semiconductors [7] . 

Compton profiles or three-dimensional electron momentum densities have been 
calculated in the past with the DFT using, for example, the full-potential 
LAPW (FLAPW) method [8,9,10,11], localized basis sets [12,13], pseudowave 
functions [14], the orthogonalized plane-wave (OPW) method [15], the re- 
construction of AE wave functions from pseudowave functions [16], the lin- 
ear muffin-tin orbital method (LMTO) [17,18,19], the full-potential LMTO 
method [20], and the Korringa-Kohn-Rostoker (KKR) band-structure scheme [21,22,23]. 
The Hartree-Fock method has also been applied [12,24], and the GW approx- 
imation [9,10,25], as well as quantum Monte Carlo (QMC) methods [26,27] 
which go beyond the DFT and the LDA. The PAW scheme has obvious ad- 
vantages compared with many of these schemes. It is formally simple and both 
efficient and fiexible in practice. First-row elements, transition metals and rare- 
earth elements can be treated. It uses the full potential and charge density, 
and AE wave functions are available. The LMTO and KKR approaches are 
usually limited to the muffin-tin form of the potential, and the pseudopotential 
method suffers from transferability problems and from the loss of information 
on the true wave functions near nuclei. The localized basis sets may have con- 
vergence problems in the interstitial regions of solids, whereas the basis set of 
the PAW method is complete and takes into account the orthogonality of the 
valence states to the core states automatically without having to treat also 
the core electrons self-consistently. Formally the closest resemblance is to the 
OPW method, which also contains a finear transformation between PS and AE 
wave functions. Also the reconstruction of AE wave functions from PS wave 
functions has some similarities, e.g. the existence of soft PS wave functions 
and the frozen-core approximation used. However, the PAW transformation is 
more effective than the OPW transformation, and the reconstruction proce- 
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dure of AE wave functions from PS wave functions is a more complicated and 
artificial way to obtain the AE wave functions than the built-in transformation 
of the PAW method. 

In this work we focus on the PAW method as a computational tool and demon- 
strate its applicability to the calculation of the valence electron momentum 
densities by calculating Compton profiles for different systems. According to 
our knowledge, this is the first time the PAW method is applied to this prob- 
lem. Ishibashi [28] has recently used the PAW method to calculate momentum 
densities of annihilating electron-positron pairs in bulk materials, and Rum- 
mukainen et al. [29] those in the case of defects is solids. We show examples 
of semiconducting, metallic and molecular systems, namely bulk Si, bulk Li, 
bulk Al, and the water dimer (H20)2, and compare our results with experiment 
and other computational results. We restrict ourselves only to the non-spin- 
polarized case, but stress that the PAW method can also be applied to the 
calculation of spin-dependent momentum densities, i.e. magnetic Compton 
profiles [30,31,32], within the impulse approximation. 



2 Method of calculation 



2. 1 Projector augmented-wave functions 



The PAW method is based on a linear transformation T between soft PS 
valence wave functions |^) and corresponding AE valence wave functions l^*), 

|*) = r|^). (1) 



The transformation differs from unity by local atom-centered contributions Tr 
such that 

r = i + ^T«. (2) 

R 



Each local contribution represents the difference between the AE wave func- 
tion and the PS wave function and acts within some augmentation region VLr 
enclosing one atom. The AE wave function \^) is expanded locally within ev- 
ery VIr into solutions of the Schrodinger equation for the corresponding 
isolated atom. The index i refers to the site index i?, the angular momentum 
indices and an additional index k referring to the reference energy e^i- 

Also the PS wave function |^) is expanded locally using a complete set of 
soft, fixed PS partial waves \4>i) and the same expansion coefficients as for 
the AE wave functions. There is exactly one PS partial wave for each AE 
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partial wave |f?)j). The choice of the PS partial waves determines the local 
contribution of to the transformation operator T, i.e., 

14) = (l+rR)|0,), within Q^. (3) 



For the transformation T to be linear, the expansion coefficients have to be 
linear functionals of the PS wave functions, i.e., inner products with 
suitable localized and fixed projector functions The transformation T can 
now be written as 

^=1 + E(I'^^)-|'^^»(P^I- (4) 

i 

The final expression for the AE wave function is then 

m-\^)+l:m-\mp^\^)■ (5) 



The total energy functional is modified using Eq. (5). The PS wave functions 
1^') are the variational quantities in the PAW method, which leads to a PS 
Hamilton operator H acting on the PS wave functions |^). The soft PS wave 
functions obtained can be represented with a modest number of plane-waves, 
because there is no norm-conservation requirement for them and the transfor- 
mation T leads to the modified orthogonality condition (^„|T'^T|\^m) = ^nm- 
The transformation T incorporates the effects of the rapid oscillations of the 
AE wave functions within the augmentation regions VIr into every expression 
in the PAW method. 

The AE partial waves have to be orthogonalized to the core states when 
necessary so that the transformation T produces only wave functions orthog- 
onal to the core electrons. The soft PS partial waves must match the AE 
partial waves outside VLr. The projector functions {pi\ are chosen so that 
the completeness condition \'Pi){Pi\ = 1 is fulfilled within every VLr. This 
enables one to expand the PS wave function l^') into PS partial waves 
locally within every VLr. The above condition implies that {pi\(i>j) = Sij. The 
construction of the partial waves and projector functions used in this work is 
described in detail in Ref. [33]. In practice, the partial wave expansions are 
truncated. Two reference energies e^i are used for each angular momentum 
quantum number / = 0,1 (/ = 0,1,2 for heavy alkali, alkali earth, and d 
elements) . 

Further details of the PAW method can be found in Refs. [4] and [34]. The 
implementation used in this work is based on the plane-wave code Vienna Ab- 
initio Simulation Package [35,36] (VASP) and its implementation of the PAW 
method [33]. 
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2.2 PAW valence electron momentum density 



The calculation of the momentum density of the AE wave functions |^) is 
straightforward. The PS wave function part |^) of Eq. (5) is by default rep- 
resented on a regular grid in the Fourier space. Note that |^) extends also to 
the augmentation regions VLr., which simplifies the calculation in comparison 
to the LAPW method. The localized partial waves are functions on a radial 
grid multiplied by spherical harmonics l^'"(f). The calculation of their Fourier 
transformations is simple, 



(|^^=l^kR{r) = 0/fc(|r - R|)i;"(r - R) 
^ cl^^=lmkR{v) = e-^P-^(/.;,(|p|) VXp)- (6) 
The radial part in the Fourier space can be written as 

1 °° 

't'l^^P) = 72^^""^"'^' / r2<^ife(r)jKpr), (7) 



where ji{pr) is a spherical Bessel function. One is now able to write down the 
expansion of the AE wavefunction into plane-waves, 

V;,k(r) = (r|*,k) = ^ E C.-k(G) exp(z(k + G) • r), (8) 

G 



where the expansion coefficient is 

C,k(G) = Qk(G) + E(</'i(k + G) - ^i(k + G))(p,|*,k)- (9) 



Here j is the band index and k the Bloch wave vector of the state, Q is the vol- 
ume of the supcrccU, G's are the reciprocal lattice points (of the superlattice) , 
and (7jk(G)'s the plane- wave expansion coefficients of the corresponding PS 
wave function |^jk)- The coefficients (pi|^jk) in Eq. (9) can be evaluated either 
in the Fourier space or in the real space depending on the implementation of 
the code used in the computation. Because the plane- waves are eigenfunctions 
of the momentum, the momentum density of the state of Eq. (8) is discrete so 
that the momentum has only values k -|- G with the probabilities determined 
by |C,k(G)p. 

In the calculation of electron momentum densities it is not required to increase 
the kinetic energy cutoff of the PS wave functions | from typical values of 
about 250 — 400 eV. The convergence of the AE wave function requires only 
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the convergence of the PS wave function. The summation in Eq. (5) converges 
when the expansion coefficients {pi\'^) have converged and the partial waves 
form a complete set of functions. Furthermore, it can be shown that the PAW 
basis set is complete whenever the plane waves form a complete set, irrespec- 
tive of the partial- wave truncation [4]. The truncation, however, affects the 
core-valence orthogonality. It also causes errors in the expectation values and 
the total energy. In Eq. (9) the contribution of the PS wave function term 
C'jk(G) to the momentum density extends up to the momentum value Gmax 
corresponding the kinetic energy cutoff of the PS wave function. The effect of 
the partial waves can be taken into account up to an arbitrary value G^^. 

We have tested the PAW basis set by calculating valence Compton profiles for 
several different isolated atoms and compared our results with the results of an 
AE code for atoms. We get a perfect match also at high momenta {pz > 5 a.u.). 
The agreement improves systematically with increasing G^^. 



2.3 Compton profiles 



In the Compton scattering experiment the quantity of interest is the so-called 
Compton profile [1] 

J{Pz)=j j dp:cdpyp{p), (10) 



where p(p) is the three-dimensional ground-state momentum density given 
within the independent-particle model by 

(11) 



^^^^ " (2^ ^ / ^'^®^P(~^P ■ 



where '0j(r)'s are the occupied single-particle states. Eq. (10) is based on 
the impulse approximation [2]. It is assumed to be valid when the energy 
transferred in the scattering process is much larger than the binding energy 
of the electronic states involved. 

In practice we calculate the momentum density of Eq. (11) by first calculat- 
ing the wavefunctions using a uniform k-point mesh dense enough to give a 
sufficient resolution. Then we obtain the momentum density by summing up 
the squares of the Fourier coefficients as 

p(p) = ^EEEIc.><(G)r4H-G,p. (12) 

" k j G 
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The standard way to take into account correlation effects beyond the independent- 
particle model in the DFT is the so-called Lam-Platzman (LP) correlation cor- 
rection [37], which is isotropic within the LDA. In general, however, compar- 
isons have shown that there are systematic anisotropic discrepancies between 
theory and experiment. As already noted, momentum densities have been cal- 
culated using the QMC methods [26,27] and the GW approximation [9,10,25]. 
However, in the case of Si [26], and Li [27] (QMC) and Li [38,25] (GW) the 
correlation correction is not found to differ appreciably from the LP corrected 
LDA results [31]. Wakoh et al. [39] have shown that the theoretical Comp- 
ton profiles of Al can be fitted to the experiment using a phenomenological 
energy-dependent occupation number instead of the LP correction. Recently, 
Barbiellini and Bansil [40] have suggested a BCS-like approach in which the 
many-body wavefunction is constructed from singlet pair wavefunctions (gem- 
inals) by taking an antisymmetrized geminal product (AGP). The scheme is 
promising in being rather successful in explaining discrepancies between exper- 
iment and previous correction schemes [31]. Corrections to both the impulse 
approximation and to the independent particle model have been recently pro- 
posed within the Dyson orbital formalism [41]. 



2.4 Computational details 



We use in our calculations the LDA based on the Ceperley-Adler [42] elec- 
tron gas results as parametrized by Perdew and Zunger [43], except in the 
case of the H2O dimer for which a gradient-corrected functional [44] is used. 
We model the structures using cubic supercells. The Brillouin zone is sam- 
pled using uniform F-centered k-point meshes. For the cubic solids we use 
the irreducible l/48th part of the Brillouin zone. We employ the frozen-core 
approximation and calculate only the valence electron momentum density. 
The core electron momentum density is isotropic and thus does not affect 
anisotropy curves. However, in the case of the bulk Li also the core electrons 
are treated self-consistently and taken into account in the Compton profiles. 
We do not include the LP correction or any anisotropic correlation correction 
in the calculated Compton profiles. Since the LP correction is isotropic within 
the LDA its effect cancels out in the anisotropy curves. The Compton pro- 
files are calculated by integrating the three-dimensional electron momentum 
density directly over planes in the p-space. The theoretical profiles are convo- 
luted with a Gaussian function with a full width at half maximum (FWHM) 
corresponding to experimental resolution. 
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3 Results and discussion 



3.1 Bulk Si 



Wc begin by studying the anisotropy of Compton profiles of Si and compare 
our results with experiment and a previous computational study. Delaney et 
al. [16] have compared the results of their PS wave function calculations to 
experimental Compton profiles. Further, they have reconstructed AE valence 
wave functions from the PS valence wave functions and studied the effect of 
the reconstruction on the Compton profiles obtained. Delaney et al. used the 
LDA, norm-conserving pseudopotentials [45] and the reconstruction scheme by 
Meyer et al. [46] . They noticed that the reconstruction improves the agreement 
with experiment substantially. The changes in the momentum density caused 
by the reconstruction are not spherically- symmetric and affect the anisotropy 
between Compton profiles along different directions improving the agreement 
with the experiment. We make the same observation and get practically iden- 
tical results with Delaney et al. when we calculate anisotropies using the AE 
and the PS wave functions of the PAW method. For the bulk Si we use the ex- 
perimental lattice constant of 10.26 a.u. The k-point mesh used is a F-centered 
10 X 10 X 10 mesh. The sampling corresponds to the momentum resolution of 
0.061 a.u. The kinetic energy cutoff for the PS wave functions is 307 eV. The 
three-dimensional momentum density used in the calculation of the Compton 
profiles [Eq. (12)] covers momenta up to Pmax = 5.0 a.u. 

Figure 1 shows our results for the [100]-[111] anisotropy compared with the 
experimental result from Ref. [47]. The transformation T from the PS wave 
functions to the AE wave functions reduces the overestimation of the am- 
plitudes of the peaks in the anisotropy. The [100]-[110] anisotropy is also 
improved but relatively less affected by the reconstruction. 



3.2 Bulk Li 



Bulk Li has been the subject of many high-resolution Compton scattering 
studies. Its Fermi surface has been studied by looking at the first and second 
derivatives of directional Compton profiles [48] and by reconstructing the full 
three-dimensional momentum density [49,50]. We compare our results for Li 
with experimental high-resolution Compton profiles and the theoretical re- 
sults from Ref. [48]. The theoretical results in Ref. [48] have been calculated 
with the KKR band structure scheme and the LDA. The momentum density 
of Li has also been calculated using the LMTO method [17,18], the FLAPW 
method [9,10,11], the GW approximation [9,10,25], and quantum Monte Carlo 
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methods [27]. Theoretical calculations based on the LDA have failed to de- 
scribe the momentum density of Li also when the Lam-Platzman correlation 
correction has been taken into account [48]. The discrepancies have been at- 
tributed to electronic correlation beyond the LDA [48,49] or to the thermal 
disorder [18] . Our calculation does not properly take into account these effects. 
Therefore, instead of experimental data we have to use other computational 
results based on the LDA as benchmarks. 

We calculate our results for the bulk Li using the experimental lattice constant 
of 6.60 a.u. The k-point mesh used is a F-centered 30 x 30 x 30 mesh. The 
sampling corresponds to the momentum resolution of 0.032 a.u. The kinetic 

energy cutoff of the PS wave functions is 271 eV. Momenta up to Pmax — 
10.0 a.u. are taken into account when calculating the Compton profiles. 

Figure 2 shows in addition to our results experimental anisotropies [48] and 
those calculated with the KKR [48]. The PAW and the KKR results are in 
perfect agreement although the mufiin-tin approximation has been used in 
the KKR calculation. A disagreement is found between our results and the 
FLAPW results by Baruah et al. [11] although in principle the only major 
difference between the calculations is the choice of the basis set. Actually, 
Baruah et al. used gradient corrections to the exchange-correlation potential, 
but concluded that their effect on the electron momentum density is not sig- 
nificant. However, the FLAPW results by Kubo [9] are in good agreement 
with the KKR method, and thus also agree with our calculation. 



3.3 Bulk Al 



Ohata et al. [51] have measured directional high- resolution Compton profiles 
of Al and compared them with theoretical results calculated with the KKR 
band structure scheme within the LDA. For Al the agreement between theory 
and experiment was found to be better than in the case of Li. Correlation 
effects are less important as the electron density is higher. Al is also harder 
than Li and other alkali metals, which reduces the effect of thermal disorder 
on Compton profiles. However, because the anisotropy of Compton profiles of 
Al is very small, and the anisotropy curves and derivatives of Compton profiles 
contain detailed features, Al is the most challenging test for the PAW method 
presented in this paper. 

For the bulk Al we use the experimental lattice constant of 7.65 a.u. The k- 

point mesh used is a F-centered 60 x 60 x 60 mesh. The sampling corresponds 
to the momentum resolution of 0.014 a.u. The kinetic energy cutoff for the 
PS wave functions is 296 eV. Momenta up to pmax = 5.0 a.u. are taken into 
account when calculating the Compton profiles. 
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Figure 3 shows the calculated valence Compton profiles and their first and sec- 
ond derivatives compared with the experimental results of Ref. [51]. As can be 
expected, for the LDA, the theoretical Compton profiles are higher at low mo- 
menta [48,16,51]. The Lam-Platzman correlation correction would shift weight 
from lower to higher momenta improving the agreement with the experiment. 
The first and second derivatives of the Compton profiles calculated with the 
PAW method agree better with experiment than the KKR results by Ohata et 
al. (see the data in Ref. [51]). Anisotropics of the Compton profiles are studied 
in Fig. 4. Agreement with the KKR scheme is not as good as in the case of 
Li. Compared with the experiment the agreement is about as good for both 
theoretical methods. The PAW method does not reproduce all the features 
of the experimental anisotropics as clearly as the KKR scheme. However, it 
predicts the positions of the peaks in the anisotropics better, especially in the 
[110]-[100] and [111]-[100] differences. 



34 Water dimer (H20)2 



Finally, we calculate the Compton profile of the gas-phase water dimer using 
the generalized-gradient approximation (GGA) by Perdew and Wang [44], 
and the experimental geometry by Dyke et al. [52] with the 0-0 separation of 
2.98 A. The water dimer can be considered as a prototype system to study how 
the constitutive interactions of the hydrogen bond (e.g. exchange repulsion, 
charge transfer and polarization) are reflected in the Compton profile [12]. The 
gas-phase dimer represents a moderately hydrogen-bonded system, where the 
main oscillatory signal in the Compton profile is expected to stem from the 
exchange interaction of the overlapping molecular wave functions [12]. Since 
the binding compared to the metallic and covalent systems is much weaker, 
the wave functions of the two water molecules conserve to a large extent their 
unperturbed form. Hydrogen bonds provide thus one of the limiting cases to 
study the applicability of the PAW calculations to describe chemical bonding. 
The PAW result is calculated using a periodic supercell with a side length of 
35.0 a.u. in order to minimize the interaction between the periodic images of 
the dimer. The k-point samphng (a F-centered 3x3x3 mesh) used when 
calculating the momentum density corresponds to the momentum resolution 
of 0.060 a.u. In the calculation of the effective potential, however, only the F 
point is used. The kinetic energy cutoff of the PS wave functions is 500 cV. 
The momentum density is taken into account up to the momentum value 
of Pmax — 13.7 a.u. when calculating the Compton profile. We calculate the 
Compton profile of the H2O dimer also by employing localized basis functions. 
The same GGA functionals are used. For oxygen a triple-zeta valence plus 
polarization type basis set is used, and for hydrogen a primitive set augmented 
by one p- function in a [3s, Ip] contraction (for more details, see Ref. [13]). No 
experimental data is available for the directional anisotropics of the dimer. 
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Figure 5 shows the anisotropy between the profile along the 0-0 direction 
(direction of the H-bond) and that in the direction orthogonal to the mirror 
plane of the H2O dimer. The solid line denotes the AE PAW result, the dashed 
line the result calculated with locahzed basis functions, and the dotted hue the 
result calculated using the PS wave functions of the PAW method. Agreement 
between the two AE calculations is very good, especially at high momenta. The 
positions of the peaks are improved discernibly when the PS wave functions 
are replaced by the AE wave functions. At momenta Pz > 2.5 a.u. the PS 
wave functions fail to describe the positions of the peaks properly. At lower 
momenta the deficiencies in the basis sets and the residual interaction between 
the periodic images in the PAW method explain the slight differences in the 
amplitudes between the AE calculations. 



4 Summciry and conclusions 

We have presented valence electron Compton profiles calculated within the 
DFT using the projector augmented- wave (PAW) method. The accuracy of the 
method has been demonstrated to be competitive with other band-structure 
and basis set schemes. The PAW method is ideal for the calculation of va- 
lence electron momentum densities of covalent systems. Also for systems with 
much weaker orbital rearrangements, such as hydrogen bonds, the PAW basis 
set yields reliable results. It becomes highly beneficial to use a unified ap- 
proach with a well-controlled basis set and full-potential treatment for the 
wide variety of nonmetallic systems. The effect of the frozen-core approxima- 
tion is negligible in most systems. The quahty of the PAW calculations can 
be systematically improved by increasing the energy cutoff parameter. For 
metallic systems (Li, Al) in the LDA, despite the coarser k-mesh, PAW yields 
essentially as good results as the KKR method. Nevertheless, for metals both 
methods suffer from discrepancies as compared with experiment, attributed 
to electron correlations beyond the LDA or to thermal disorder. Our work 
suggests that the PAW method which is becoming ever more popular in elec- 
tronic structure calculations suits well for describing a wide variety of ordered 
and disordered covalently or hydrogen bonded systems. 
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Fig. 1. Bulk Si. The anisotropy between the Compton profiles along the [100] and 
[111] directions. The theoretical curves arc convoluted with a Gaussian function 
with a FWHM of 0.13 a.u. corresponding to the experimental resolution. The ex- 
perimental data is from Ref. [47] . 
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Fig. 2. Bulk Li. Anisotropics between Compton profiles along different directions. 
The theoretical curves are convoluted with a Gaussian function with a FWHM of 
0.12 a.u. corresponding to the experimental resolution. The experimental data and 
the KKR result are from Ref. [48]. 
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Fig. 3. Valence Compton profiles of Al along different directions (the first row). Also 
their first and second derivatives obtained by numerical differentiation are shown 
(the second and third rows). The theoretical Compton profiles are convoluted with a 
FWHM of 0.12 a.u. corresponding to the experimental resolution. The experimental 
valence Compton profiles (from Ref. [51]) have been obtained by subtracting the 
KKR core profile from the measured ones. 
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Fig. 4. Bulk Al. Anisotropies between Compton profiles along different directions. 
The theoretical curves are convoluted with a Gaussian function with a FWHM of 
0.12 a.u. corresponding to the experimental resolution. The experimental data is 
from Ref. [51]. 
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